R jako GIS
Podpůrný text pro předmět: GIS pro biologické aplikace
Autor: Matěj Man
Aktualizace: 18. 11. 2019
Odkaz na R skript: Kopírovat a vložit do R
Další zdroje:
- Online kniha Geocomputation with R
- Online kniha Spatial Data Science with R
- Velice pěkný návod pro R knihovnu sf
- The Visual Raster Cheat Sheet Rpubs dokument
Knihovny
list.of.packages <- c("sf","raster","mapview","whitebox","randomcoloR","leaflet")
# install.packages(list.of.packages)
library(sf)
library(raster)
library(mapview)
library(randomcoloR)
library(leaflet)Načítání vektorových GIS dat do R
## Načítání vektorových dat shp
# nastavte kde máte u sebe na PC data
# pozor musíte zdvojit nebo otočit lomítka
cesta<-"d:\\Owncloud\\ŠKOLA\\Učení\\GIS\\cv 6\\2020_2021_cv6\\"
# zkonstuuje cestu kde leží vrstva Brdy
data.path<-paste0(cesta,"\\parky2\\CHKO_Brdy.shp")
# načte .shp do R
brdy<-st_read(data.path,stringsAsFactors = F,quiet = T)
# zobrazí načtený soubor v interkativním okně mapy
# mapview(brdy)
# obdobně načteme třeba hranici ČR
data.path<-paste0(cesta,"\\hrcr1_wgs.shp")
hrcr<-st_read(data.path,stringsAsFactors = F,quiet = T)
# prostý obrázek, bez interaktivity
plot(st_geometry(hrcr))
# parametr add přidá další vrstvu do existujícího obrázku
plot(st_geometry(brdy),col="red",add=T) 
# Nastaví pracovní adresář working dorectory (WD)
# dále nemusíme data z WD volat celou cestou, stačí názvy
setwd(cesta)
## Načítání vektorových dat z tabulky
# načíst prostou tabulku
chmu<-read.table("staniceCHMUtablecoma.csv",header = T, sep=",")
# převést tabulky na prostorová data
chmu_sf<-st_as_sf(chmu, coords = c("Xcoo", "Ycoo"),crs = 4326)
## vizualizovat prostorová data interaktivně
# pokročilejší balíček leaflet
# definujeme paletu o 17 barvách podle typu stanice
distCol<- colorFactor(distinctColorPalette(17), chmu_sf$Station.ty)
# definoujeme okno s mapou
leaflet(chmu_sf) %>%
addTiles() %>%
addCircleMarkers(popup=chmu_sf$Name,color = ~distCol(Station.ty),
stroke = FALSE, fillOpacity = 0.8,radius=4) %>%
addLegend(pal = distCol, values = ~Station.ty)Projekce a crs
st_crs(brdy) # jaký crs má vrstva brdy ?## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]st_crs(hrcr) # jaký crs má vrstva hrnice čr ?## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]## vektory
# transformace podle EPSG
brdy_32633<-st_transform(brdy,32633)
st_crs(brdy_32633) # jaký crs má vrstva brdy_32633 ?## Coordinate Reference System:
## User input: EPSG:32633
## wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 33N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["World - N hemisphere - 12°E to 18°E - by country"],
## BBOX[0,12,84,18]],
## ID["EPSG",32633]]# transformace podle proj4 string
brdy_5514<-st_transform(brdy, "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
st_crs(brdy_5514) # jaký crs má vrstva brdy_5514 ?## Coordinate Reference System:
## User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
## wkt:
## BOUNDCRS[
## SOURCECRS[
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["Unknown based on Bessel 1841 ellipsoid",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]],
## CONVERSION["unknown",
## METHOD["Krovak (North Orientated)",
## ID["EPSG",1041]],
## PARAMETER["Latitude of projection centre",49.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8811]],
## PARAMETER["Longitude of origin",24.8333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8833]],
## PARAMETER["Co-latitude of cone axis",30.2881397222222,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",1036]],
## PARAMETER["Latitude of pseudo standard parallel",78.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8818]],
## PARAMETER["Scale factor on pseudo standard parallel",0.9999,
## SCALEUNIT["unity",1],
## ID["EPSG",8819]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Position Vector transformation (geog2D domain)",
## ID["EPSG",9606]],
## PARAMETER["X-axis translation",570.8,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",85.7,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",462.8,
## ID["EPSG",8607]],
## PARAMETER["X-axis rotation",4.998,
## ID["EPSG",8608]],
## PARAMETER["Y-axis rotation",1.587,
## ID["EPSG",8609]],
## PARAMETER["Z-axis rotation",5.261,
## ID["EPSG",8610]],
## PARAMETER["Scale difference",1.00000356,
## ID["EPSG",8611]]]]# transformace podle jiné existující vrstvy
brdy_krovak<-st_transform(brdy, st_crs(brdy_5514))
st_crs(brdy_krovak) # jaký crs má vrstva brdy_5514 ?## Coordinate Reference System:
## User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
## wkt:
## BOUNDCRS[
## SOURCECRS[
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["Unknown based on Bessel 1841 ellipsoid",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]],
## CONVERSION["unknown",
## METHOD["Krovak (North Orientated)",
## ID["EPSG",1041]],
## PARAMETER["Latitude of projection centre",49.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8811]],
## PARAMETER["Longitude of origin",24.8333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8833]],
## PARAMETER["Co-latitude of cone axis",30.2881397222222,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",1036]],
## PARAMETER["Latitude of pseudo standard parallel",78.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8818]],
## PARAMETER["Scale factor on pseudo standard parallel",0.9999,
## SCALEUNIT["unity",1],
## ID["EPSG",8819]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Position Vector transformation (geog2D domain)",
## ID["EPSG",9606]],
## PARAMETER["X-axis translation",570.8,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",85.7,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",462.8,
## ID["EPSG",8607]],
## PARAMETER["X-axis rotation",4.998,
## ID["EPSG",8608]],
## PARAMETER["Y-axis rotation",1.587,
## ID["EPSG",8609]],
## PARAMETER["Z-axis rotation",5.261,
## ID["EPSG",8610]],
## PARAMETER["Scale difference",1.00000356,
## ID["EPSG",8611]]]]# kde jsou ty brdy? no nevidím je protože to má jiný CRS
plot(st_geometry(hrcr),main="Kde jsou ty Brdy?")
plot(st_geometry(brdy_32633),add=T) 
# musíme tedy transformovat jedno nebo druhé
hrcr_32633<-st_transform(hrcr, 32633)
plot(st_geometry(hrcr_32633),main="No jo, už máme správný CRS")
plot(st_geometry(brdy_32633),add=T,col="red")
Načítání rastrových GIS dat do R
# načítání jdnoduchou funkcí "raster"
DEM<-raster(paste0(cesta,"\\DEM_Jested_100m.tif"))
# kontrolí obrázek
# data tu jsou
plot(DEM)
# ale leží na správném místě?
# ukáže crs vrstvy
DEM@crs## CRS arguments:
## +proj=krovak +axis=swu +lat_0=49.5 +lon_0=24.8333333333333 +alpha=0
## +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs# zkušeným okem, nebo podle errorů při další alanýze poznáte,
# že nemá dobře definovaný crs, je to křovák, ale šptně definovaný.
# tak mu řekneme, jaká je korektní definice křováka
# http://freegis.fsv.cvut.cz/gwiki/S-JTSK
DEM@crs<-crs("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
# i rastr je možné vizualizovat interaktivně
# definujeme barvy
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(DEM),
na.color = "transparent")
leaflet(DEM) %>%
addTiles() %>%
addRasterImage(DEM,colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(DEM),
title = "Elevation [m]")Příklad vektorové analýzy
# plocha polygonů
st_area(brdy)## 345041249 [m^2]# délka linií (obvod polygonů)
st_length(brdy)## 195793.8 [m]# buffer 5 km
brdy_5000<-st_buffer(brdy_32633,5000)
plot(st_geometry(brdy_5000), main="buffer 5 km")
plot(st_geometry(brdy_32633),add=T,col="red")
# centroidy
brdy_cen<-st_centroid(brdy_32633) ## Warning in st_centroid.sf(brdy_32633): st_centroid assumes attributes are
## constant over geometries of xplot(st_geometry(brdy_32633), main="centroid",graticule=T,axes=T)
plot(st_geometry(brdy_cen),add=T,col="blue",pch=19)
# průnik
chmu_brdy<-st_intersection(brdy,chmu_sf)## although coordinates are longitude/latitude, st_intersection assumes that they are planar## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(st_geometry(brdy), main="intersection - chmu stanice v brdech",graticule=T,axes=T)
plot(st_geometry(chmu_brdy),add=T,col="red")
Hromadné zpracování dat
Aneb co by nám v QGIS trvalo klikat hodiny můžeme v R naprogramovat během minut
# načteme všechny cesty k shp vrstvám v daném adresáři
parky.files<-list.files(path=paste0(cesta,"\\parky2"),pattern="*.shp$",full.names = T)
## když chceme všechny soubor do jednoho objektu
# načte první shpfile
parky.init<-st_read(parky.files[1],quiet = TRUE)
# připojí všechny další shp file do jendoho objektu
for (i in 2:length(parky.files)) {
parky.next<-st_read(parky.files[i],quiet = TRUE)
parky.init<-rbind(parky.init,parky.next)
}
# připojit metadat GIS join
meta.path<-paste0(cesta,"\\metadata_parky.csv")
metadata<-read.table(meta.path,sep=";",header=T,stringsAsFactors = T)
parky<-merge(parky.init,metadata)
# kontrolní obrázek
plot(st_geometry(hrcr))
plot(st_geometry(parky),add=T,col="green")
# interaktivně
distCol<- colorFactor(distinctColorPalette(32), parky.init$OBJECTID)
leaflet(parky) %>%
addTiles() %>%
addPolygons(color = ~distCol(OBJECTID),
stroke = FALSE, fillOpacity = 0.8,
label = paste0(parky$KAT,"_",parky$NAZEV))Rastrové analýzy
Velice dobrá knihovna pro rastrové analýzy v R je whitebox, což je implementace jinak v Python/Rust napsaných geo algoritmů.
## instalace balíčku whitebox pro rastrové analýzy
# install.packages("whitebox", repos="http://R-Forge.R-project.org") # instalace
# whitebox::wbt_init() # inicializace
library(whitebox)
print(wbt_help())## [1] "WhiteboxTools Help"
## [2] ""
## [3] "The following commands are recognized:"
## [4] "--cd, --wd Changes the working directory; used in conjunction with --run flag."
## [5] "-h, --help Prints help information."
## [6] "-l, --license Prints the whitebox-tools license."
## [7] "--listtools Lists all available tools. Keywords may also be used, --listtools slope."
## [8] "-r, --run Runs a tool; used in conjuction with --wd flag; -r=\"LidarInfo\"."
## [9] "--toolbox Prints the toolbox associated with a tool; --toolbox=Slope."
## [10] "--toolhelp Prints the help associated with a tool; --toolhelp=\"LidarInfo\"."
## [11] "--toolparameters Prints the parameters (in json form) for a specific tool; --toolparameters=\"LidarInfo\"."
## [12] "-v Verbose mode. Without this flag, tool outputs will not be printed."
## [13] "--viewcode Opens the source code of a tool in a web browser; --viewcode=\"LidarInfo\"."
## [14] "--version Prints the version information."
## [15] ""
## [16] "Example Usage:"
## [17] ">> .\\whitebox-tools.exe -r=lidar_info --cd=\"\\path\\to\\data\\\" -i=input.las --vlr --geokeys"
## [18] ""# pozor na diakritiku a mezery v cestě. Pro WBT nesmí být
# cesta k rastrovým datům
dem <-("c:\\Users\\manma\\Downloads\\DEM_Jested_100m.tif")
## stínovaný reliéf
# cesta k budoucímu výsledku
output<- ("c:\\Users\\manma\\Downloads\\Hillshade_100m.tif")
# spustí algoitmus pro výpočet stinovaného reliéfu
wbt_hillshade(dem, output, azimuth = 315, altitude = 30,
zfactor = 1,verbose_mode = FALSE)## [1] "hillshade - Elapsed Time (excluding I/O): 0.4s"# načtu výsledek jako rastr
hill<-raster(output)
# načtu původní dem jako rastr
elev<-raster(dem)
# vytisknu oba obrázky pro kontrolu
plot(hill, col=grey(0:100/100), legend=FALSE,main="stínovaný reliéf")
plot(elev, col=rainbow(25, alpha=0.35), add=TRUE)
## Sklon svahu
output<- ("c:\\Users\\manma\\Downloads\\Slope_100m.tif")
wbt_slope(dem,output,zfactor = 1, verbose_mode = FALSE)## [1] "slope - Elapsed Time (excluding I/O): 0.3s"slp<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="Sklon svahu")
plot(slp, col=heat.colors(25, alpha=0.2), add=TRUE)
## The terrain ruggedness index (TRI)
output<- ("c:\\Users\\manma\\Downloads\\TRI_100m.tif")
wbt_ruggedness_index(dem,output,zfactor = 1, verbose_mode = FALSE)## [1] "ruggedness_index - Elapsed Time (excluding I/O): 0.4s"tri<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="terrain ruggedness index")
plot(tri, col=cm.colors(25, alpha=0.3), add=TRUE)